#include <stdio.h>
#include <stdlib.h>

#include <math.h>

int fs_sexa (char *out, double a, int w, int fracbase)
{
	char *out0 = out;
	unsigned long n;
	int d;
	int f;
	int m;
	int s;
	int isneg;

	/* save whether it's negative but do all the rest with a positive */
	isneg = (a < 0);
	if (isneg)
	    a = -a;

	/* convert to an integral number of whole portions */
	n = (unsigned long)(a * fracbase + 0.5);
	d = n/fracbase;
	f = n%fracbase;
	/* form the whole part; "negative 0" is a special case */
	if (isneg && d == 0)
	    out += sprintf (out, "%*s-0", w-2, "");
	else
	    out += sprintf (out, "%*d", w, isneg ? -d : d);

	/* do the rest */
	switch (fracbase) {
	case 60:	/* dd:mm */
	    m = f/(fracbase/60);
	    out += sprintf (out, ":%02d", m);
	    break;
	case 600:	/* dd:mm.m */
	    out += sprintf (out, ":%02d.%1d", f/10, f%10);
	    break;
	case 3600:	/* dd:mm:ss */
	    m = f/(fracbase/60);
	    s = f%(fracbase/60);
	    out += sprintf (out, ":%02d:%02d", m, s);
	    break;
	case 36000:	/* dd:mm:ss.s*/
	    m = f/(fracbase/60);
	    s = f%(fracbase/60);
	    out += sprintf (out, ":%02d:%02d.%1d", m, s/10, s%10);
	    break;
	case 360000:	/* dd:mm:ss.ss */
	    m = f/(fracbase/60);
	    s = f%(fracbase/60);
	    out += sprintf (out, ":%02d:%02d.%02d", m, s/100, s%100);
	    break;
	default:
	    printf ("fs_sexa: unknown fracbase: %d\n", fracbase);
            return -1;
	}

	return (out - out0);
}

#define MATRIX_LOG(name, in) \
  printf("Matrix %s:\n%g %g %g\n%g %g %g\n%g %g %g\n", name, in[0][0], in[0][1], in[0][2], in[1][0], in[1][1], in[1][2], in[2][0], in[2][1], in[2][2])

void mult_matrix_3x3(double in1[3][3], double in2[3][3], double out[3][3]) {
  int i, j, k;
  for (i=0; i < 3; i++) {
    for (j=0; j < 3; j++) {
      out[i][j] = 0.0;
      for (k=0; k < 3; k++)
	out[i][j] += in1[i][k] * in2[k][j];
    }
  }
}

/* Taki' paper p43 */
#define DEFAULT_THETA 4.0
#define DEFAULT_GAMMA 3.0
#define DEFAULT_LATITUDE 49.4944

int main(int argc, char** argv)
{
  double theta = DEFAULT_THETA * M_PI / 180.0;
  double gamma = DEFAULT_GAMMA * M_PI / 180.0;
  double tazimuth, taltitude;

  double lat = DEFAULT_LATITUDE * M_PI / 180.0;
  double colat;
  double ratel, Htel, dectel;
  double aztel, alttel;
  char convalt[13], convaz[13];

  double toaltaz[3][3];
  double Lt, Mt, Nt;
  double Lht, Mht, Nht;

  double rotz[3][3];
  double rotx[3][3];
  double mat[3][3];

  double lst=0.0;
  double ra, dec, H; 
  double Lc, Mc, Nc;

  double mra, mdec;
  double L, M, N;
  //printf("String ra %s dec %s\n", argv[1], argv[2]);
  sscanf(argv[1], "%lf", &lst);
  sscanf(argv[2], "%lf", &ra);
  sscanf(argv[3], "%lf", &dec);

  ratel = (lst + 6.0) * M_PI / 12.0;
  ratel -= theta;
  Htel = (lst * M_PI / 12.0) - ratel;
  dectel = (M_PI / 2.0) - gamma;
  colat = lat - (M_PI / 2.0);
  toaltaz[0][0]=cos(colat); toaltaz[0][1]=0.0; toaltaz[0][2]=sin(colat);
  toaltaz[1][0]=0.0; toaltaz[1][1]=1.0; toaltaz[1][2]=0.0;
  toaltaz[2][0]=-sin(colat); toaltaz[2][1]=0.0; toaltaz[2][2]=cos(colat);
  Lt = cos(dectel) * cos(-Htel);
  Mt = cos(dectel) * sin(-Htel);
  Nt = sin(dectel);
  Lht=toaltaz[0][0] * Lt + toaltaz[0][1] * Mt + toaltaz[0][2] * Nt;
  Mht=toaltaz[1][0] * Lt + toaltaz[1][1] * Mt + toaltaz[1][2] * Nt;
  Nht=toaltaz[2][0] * Lt + toaltaz[2][1] * Mt + toaltaz[2][2] * Nt;
  
  aztel = atan(Mht/Lht) * 180.0 / M_PI;
  if (Lht < 0) aztel += 180.0;
  aztel = -180.0 - aztel;
  while (aztel<0.0) aztel+=360.0;
  while (aztel>360.0) aztel-=360.0;

  alttel = asin(Nht) * 180.0 / M_PI;

  fs_sexa(convalt, alttel, 3, 3600);
  fs_sexa(convaz, aztel, 3, 3600);
  //printf("Telescope: Z drift is %g rad, X drift is %g rad\n",  theta, gamma);
  printf("Telescope: alt %s (%g) az %s (%g)\n", convalt, alttel, convaz, aztel);

  //printf("Read ra %g dec %g\n", ra, dec);
  H=(lst - ra) * M_PI / 12.0;
  dec=dec * M_PI / 180.0; 

  // Taki's error ?: theta and gamma clockwise -> indirect
  //  rotation matrix are for direct angles
  //theta = -theta;
  //gamma = -gamma;
  // No this concerns the coordinate systme where angle are reversed
  // So that's ok
  rotz[0][0]=cos(theta); rotz[0][1]=-sin(theta); rotz[0][2]=0.0;
  rotz[1][0]=sin(theta); rotz[1][1]=cos(theta); rotz[1][2]=0.0;
  rotz[2][0]=0.0; rotz[2][1]=0.0; rotz[2][2]=1.0;

  rotx[0][0]=1.0; rotx[0][1]=0.0; rotx[0][2]=0.0;
  rotx[1][0]=0.0; rotx[1][1]=cos(gamma); rotx[1][2]=-sin(gamma); 
  rotx[2][0]=0.0; rotx[2][1]=sin(gamma); rotx[2][2]=cos(gamma);

  mult_matrix_3x3(rotx, rotz, mat);
  //MATRIX_LOG("Rot x", rotx);
  //MATRIX_LOG("Rot z", rotz);
  //MATRIX_LOG("Matrix", mat);

  Lc=cos(dec) * cos(-H);
  Mc=cos(dec) * sin(-H);
  Nc=sin(dec);

  L=mat[0][0] * Lc + mat[0][1] * Mc + mat[0][2] * Nc;
  M=mat[1][0] * Lc + mat[1][1] * Mc + mat[1][2] * Nc;
  N=mat[2][0] * Lc + mat[2][1] * Mc + mat[2][2] * Nc;
  
  printf("RA %g H %g DEC %g\n", ra, H, dec);
  printf("Lc=%g Mc=%g Nc=%g\n", Lc, Mc, Nc);

  mra=atan(M/L) * 12.0 / M_PI;
  printf("atan(M/L) %g L=%g M=%g N=%g\n", mra, L, M, N);
  if (L < 0.0) mra = 12.0 + mra;
  mra+=lst;
  while (mra<0.0) mra+=24.0;
  while (mra>24.0) mra-=24.0;
  mdec=asin(N) * 180.0 / M_PI;


  printf("%.8g %.8g\n", mra, mdec);

  exit(0);
}
